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The eigendecomposition of the coupling matrix of large biological networks is central to the study 
of the dynamics of these networks. For neural networks, this matrix should reflect the topology of 
the network and conform with Dale’s law which states that a neuron can have only all excitatory or 
only all inhibitory output connections, i.e., coefficients of one column of the coupling matrix must 
all have the same sign. The eigenspectrum density has been determined before for dense matrices 
Jij [1], when several populations are considered [2, 3]. However, the expressions were derived under 
the assumption of dense connectivity, whereas neural circuits have sparse connections. Here, we 
followed mean-field approaches [4] in order to come up with exact self-consistent expressions for the 
spectrum density in the limit of sparse matrices for both symmetric and neural network matrices. 
Furthermore we introduced approximations that allow for good numerical evaluation of the density. 
Finally, we studied the phenomenology of localization properties of the eigenvectors. 


The dynamics of diverse biological systems, such as 
neural, ecological or genetic networks involves an inter¬ 
play between many individual elements. Since the precise 
nature of this coupling is difficult to determine, it is of¬ 
ten useful to consider random coupling. Specifically, in 
neuroscience synaptic connections between neurons un¬ 
derlie the dynamics of these networks, yet despite ma¬ 
jor efforts [5], the pattern of these connections is largely 
unknown. Accordingly, models of these dynamics are of¬ 
ten studied under the assumption of random connectiv¬ 
ity. The stability analysis of random networks occupies 
a central role in the study of the dynamical behavior of 
many classes of neural networks [6]. Studying the case of 
random connectivity is of further importance since it will 
serve as an important baseline for deciphering the effects 
of more specific connectivity patterns, such as structural 
motifs [7, 8]. 

Like most biological systems, realistic neuronal net¬ 
works do not have all-to-all connectivity. Instead, con¬ 
nectivity is typically highly sparse, i.e., most of the co¬ 
efficients of the connectivity matrix are zero. The 
spectrum density of Jij has been previously determined 
for dense, [1-3, 9]. Here we study the case of highly 
sparse matrices where the number of non-zero elements 
per column is finite. We are able to derive expressions 
for the eigenspectrum of sparse networks obeying the 
central demarcating line in terms of connectivity struc¬ 
ture in neural circuits, known as Dale’s law, that states 
that neural circuits are split into two populations: ex¬ 
citatory neurons whose activity evokes activity in their 
downstream neurons, and inhibitory neurons, whose ac¬ 
tivity suppresses activity in down stream neurons. We 
find striking differences both with the sparse symmetric 
case where a non-finite support tail is observed [10, 11] 
and with the dense non-symmetric case respecting Dale’s 
law in the bulk of the spectrum [1]. 

The spectrum of large, sparse, but symmetric random 
matrices has previously been studied by several meth¬ 
ods [10]. However, in most biological systems the cou¬ 


pling between units is non-symmetric. Here, we develop 
an approximation scheme based on the cavity method 
and apply it to studying the eigenvalue spectrum of non- 
symmetric matrices (the application of the method to 
symmetric matrices is outlined in [12]) and formulated 
in a different way by [13]. Next we extend our results to 
networks whose structure is non-uniform, and depends on 
the functional class of a unit, since this is the typical case 
in biological networks. Specifically, in neural circuits, the 
central demarcating line in terms of connectivity is that 
known as Dale’s law, splitting neurons into two popu¬ 
lations: excitatory neurons whose activity evokes activ¬ 
ity in their downstream neurons, and inhibitory neurons, 
whose activity suppresses activity in down stream neu¬ 
rons. Accordingly, we develop our methods below to be 
able to deal with networks composed of two neural pop¬ 
ulations. 

For general non-symmetric Jij, we follow the field- 
theoretic mapping of the eigen-spectrum density of [14, 
15] and have: 

= z*/-At) W 


where A is the Schur decomposition oi J, z = x + iy and 
9 = (5x - idy)l2, d* = (a^; - idy)l2. 

With the gaussian integral representation of the deter¬ 
minant and the proper change of basis, it follows (ignor¬ 
ing irrelevant prefactors within the log) that: 
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and integration is over the 2iV-dimensional complex field 
Ip. Note that k and k' make the integral well defined but 
are also introduced for latter convenience. The cavity 
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method will consist in isolating a small number of fields 
from the partition function Z = J d'tp evi'p {—'H). Let us 
first isolate the fields corresponding to a neuron k: 


Taking the average over the quenched variables, we 
obtain the relationship (which is asymptotically exact in 
the sparse limit): 


Z = Kk J dtpkdtpN+kdhkdhN+k exp {-T-Lk) (4) 

with: 

'Hk = K IV'fcP + k' IV-w+feP + 2m{z'iljN+kip*k) 

- 2i^{hu^Pl) - 2i^{-4,N+kh%+k) + 1/2hISfc (5) 
where we define the fields acting on and ipN+k as: 

V JkifpN+i, hN+k = Jiki’i 

( 6 ) 

Note that the dependence of Hk in hk and /lAr+fc is 
quadratic since Z is a gaussian integral and that Sfc re¬ 
mains to be determined. The variance of the external 
fields are expressed as: 

{\hk\ )\k ~2,<Thf. = ^ ' JklJkl'{'4’N+l'tpN+l')\k (7a) 
1 , 1 ' 

{\hN+k\ )\k = ^ ' JlkJl'k{'4’l'4’l')\k (7b) 

1,1' 

{hkh*j^j^^)\k = ^ JklJl'k{'ipN+li’l>)\k (7c) 

1,1' ^k 

Here the moments are computed with a hamiltonian 
lacking the terms containing 'ipk and ipN+k- The cav¬ 
ity approximation consists in neglecting the off-diagonal 
terms of these expressions: 

~'^\Jki\'^i\^N+i\'^)\k (8a) 

l^k 

~ {hkhN+k)\k « 0 (8b) 

l^k 

This approximation is thus valid when the network is 
sparsely connected (different fields are weakly correlated) 
or when the mean of the coupling constants is zero. It 
is possible following the same scheme but isolating two 
fields in the hamiltonian, to go beyond this approxima¬ 
tion and capture finite size effects at the order 1/A^. Per¬ 
forming this integral, we obtain: 

Z = (|z|2 + {k + 2alJ{K' + 

(9) 

This gives the moments: 

(l^fcP) = 2<,,/(kP+4«^J (10a) 

= 2alJ{\z\^+dalal^J (10b) 

Before deriving the eigenspectrum we extend the treat¬ 
ment, dealing with two population coupling matrices. 


{{Wk\^))) = 

^'^kl{,Jkl)d'J^lk{^Jik) 

\Jki\'^{\tpN+i\‘^)\k 

( 11 ) 

A similar expression holds for (( {\tpN+k \'^))) 

We now consider that we have two populations of neu¬ 
rons: fN columns and (1 — f)N columns of the matrix 
have coefficients independently distributed with distribu¬ 
tion TTexiJ) and 7ri„h(J) respectively. For -ipi correspond¬ 
ing to an inhibitory neuron i, the measure in the previous 
integral can now be written: 

fNNN 

dlTiYii^ {Jkl) dir^xidkl ) dTTinh (d/fc) (12) 

1=1 l=fN+l 1=1 

We have a similar expression for the ipey corresponding 
to an excitatory neuron. 

From the general expression (2), we can derive the 
spectrum density: 

Piz) = :^9^{ipk^*f^+k) (13a) 

k 

"^^?kP + 4<(kPX,,(N2) 

This expression gives the correct result in the dense limit 
as well, as shown below. 

In the dense limit, the central limit theorem holds for 
the expression (8). It is thus expected that the dense 
limit expression holds as soon as the number of connec¬ 
tion per neuron tends to infinity in the thermodynamic 
limit {ie the central limit theorem holds). In the dense 
limit, the moments are self-averaging and we can derive 
the external field variances, which for an inhibitory col¬ 
umn i gives: 

2al ^ Nifalm^+el^))) + (I - /)a?„h(((IV'^+.H))) 

( 14 ) 

We obtain: 
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and the self-consistency relationship (10) gives: 


(15a) 

(15b) 


fr(|z|2) = 0, 

|2r(|z|2) = {l + a)\z\^+ VA if jzp < z'^ 

(16) 

where A = (A^/crfnh “ l•2^P(l + + 4q;|zP(z/ - \z\'^), 

'k = and al,., af„^, are the variances of TTgx, 7ri„h 

respectively, z^ = + (1 - /)o-/„h)- 
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FIG. 1. Spectrum density of large sparse random matrices. 
The numerical evaluation comes from a direct eigendecompo- 
sition of 4000 matrices of size 10000x10000. The parameters 
are as follows: /iex = fJ-inh = 8, Jex = 1, Jinh = 1.5 and 
/ = 0.6. The dense evaluation comes from expression 17 
whereas the sparse evaluation comes from expression 20. 


We finally obtain, by plugging this result into equa¬ 
tion (13): 


. ^ ^ , (l-/)(l-a)(r(|zn-NV(|^n) 

(|zP+r)2 

(17) 

It is straightforward to verify that this expression leads 
to the same density as derived previously by Rajan and 
Abbott [1] for the dense case. 

In the sparse limit, it is no longer possible to apply the 
central limit theorem. Our approximation entails keep¬ 
ing the previous expression even in the case of sparse ma¬ 
trices whereby the self-averaging property is not strictly 
fulfilled. We begin by considering that every excitatory 
(resp. inhibitory) neuron has the same strength of cou¬ 
pling with their postsynaptic partners Je (resp. Ji) and 
the noise in the matrix only comes from the presence or 
absence of a synaptic connection: 


7rex(J') = (1 - - Jex) (18a) 

^inh(J) = (1 - /ii„h/W)5(J) + fiinb/NS{J + Jinh) 

(18b) 


We work in the sparse limit where each neuron has a 
Poissonian number of excitatory and inhibitory input 
neurons. This leads to the new expression of the self- 
consistent relationships for inhibitory neurons i and ex¬ 
citatory neurons e: 




+ 00 
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k,l,m,n —0 


p{k, /, m, n) 

NP + J?nh9ikJ)him,n) 


(19) 


where: 

p{k, l,m,n)= (0P(i-/)/x,„h im)pff,,^ (n) 

g{k,i) = k{mf)))+im,f))) 

h{m,n) = mJl^M\i^N+i\^))) + riJ^^{{{\^N+e\^))) 

We have a similar expression for (((I'i/l’eP)))- In this ex¬ 
pression, Pfj,{n) = pPe~^/n\ represents the mass of the 
Poisson distribution of parameter p. In order to apply 
the self-consistent relationship in the sparse case, we con¬ 
sidered that the neurons connected with the neuron that 
we isolated in the mean-field approximation are all equiv¬ 
alent. This is only true when the neurons are almost all 
interconnected and form a large cluster. In other words, 
the network needs to have reached the percolation three- 
shold. We show in the supplementary information [12] 
that this happens for the two-population networks when 
the average number of output connection of a neuron is 
higher than one. 

It is difficult to solve the self-consistent (19) equations 
explicitly. However, it is still straightforward to solve 
them numerically. We then obtain thanks to (13): 


P{z) 


^ V \z\^ + JL9(k,l)h{m,n) 

fq{k,l,m,n) \ 
+ JL9ik,l)h{m,n) J 


( 20 ) 


where: 

q{k,l, TO, n) = (k) pj^^^ (1) (to) pj^^^ (n) 

The comparison between this expression, the expected 
density from the dense expression, and the density ob¬ 
tained numerically is presented in figure 1. 

The sparse limit asymptotic correction to the dense 
expression are derived in [12] close to the critical value of 
\z\'^ when the spectrum density becomes zero. We show 
that the correction at that point are of order 1 /p where p 
is the mean number of outgoing connections of a neuron. 
We expect that the order of the correction is the same for 
the whole bulk of the spectrum, ie outside the divergence 
in =0. 

For comparison, we have plotted the spectrum density 
for symmetric matrices 2. It is known that for symmetric 
sparse matrices, the tail of the spectrum is finite contrary 
to the dense case. This is well captured by our approx¬ 
imation by contrast with the EM A and SDA. For non- 
symmetric matrices, the situation is different and both 
within our approximation and by numerical evaluation, 
we found that the tail of the spectrum tends to zero with 
N. 

It has been shown that for symmetric matrices, the 
non-finite spectrum support in the sparse limit is associ¬ 
ated with localized eigenvectors on the boundary of the 










4 



FIG. 2. Spectrum density of a sparse symmetric matrix. The 
different curves compare the EMA, SDA and the reported ap¬ 
proximation on the tail of the spectrum. The average number 
of connection per neuron is 10. Inset: tail of the spectrum. 


spectrum. In the non-symmetric case however, there ex¬ 
ists a critical value for the boundary of the spectrum 
identical to the one derived from the dense expression. 
We computed the participation ratio of the eigenvector 
associated with eigenvalues on the whole spectrum. This 
is displayed on figure 3. For comparison we plotted the 
same graph for the dense case in the supplementary in¬ 
formation [12]. We observe that the localization is indeed 
increased on the boundary of the spectrum, contrary to 
the dense case. However, the localization is weak and 
tends to zero with the size of the system (contrary to the 
symmetric case). 

It has been proposed in several works [16] that neural 
networks stand at the edge of chaos, close to the critical 
value derived here. We have shown that the more realis¬ 
tic assumption of sparse connections modifies the struc¬ 
ture of eigenvalues and eigenvectors around the transi¬ 
tion, and hence the states that appear close to the transi¬ 
tion. Moreover, the influence of sparsity on the tail of the 
spectrum (which control the stability behavior of the sys¬ 
tem) is very different for symmetric and non-symmetric 
matrices. This is due to the lack of feedback on more 
connected nodes for non-symmetric matrices. Our work 
hence suggests that structural motifs within neural net¬ 
works (which include over-represented reciprocal connec¬ 
tions) would have a qualitative influence on the tail of 
the eigenvalue spectrum of the coupling matrix of these 
networks. This paves the way to studying the non-linear 
dynamical behavior close to the transition in the sparse 
regime, and understanding the relationship between cir¬ 
cuit motif structure and dynamics. 

We would like to thank L. Abbott, D. Hansel, S. 
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K(2) 

FIG. 3. Eigenvalue distribution and the complex plane and 
inverse participation ratio (IPR) of the associated eigenvec¬ 
tors which represent the level of sparsity of the eigenvector 
coefficients associated to the displayed eigenvalue. The cor¬ 
responding matrix has the same parameters as the one used 
for figure 1. 


Romani and V. Samalam for useful comments on the 
manuscript. 
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SYMMETRIC MATRICES 


We consider a random symmetric matrix J of dimension N whose coefficients Jij, 1 < 
i ^ j N are identically and independently distributed with some distribution V of mean 
0 and of standard deviation a. 

Its spectrum is real with eigenvalues A*, 1 < i < iV and by dehnition its spectrum density 
as a function of /i G R is: 

1 ^ 

P(h) = Ai) (1) 

i=l 

By using the following representation of the Dirac delta distribution: 


5ix) = — limTy- 

TT e-5>0 X — le 

we obtain the following expression for the spectrum density: 

- N 




le] 


i=l 


p{^A = 


Nn 


N 

lim log / d(f) ( 

e —>-0 ^ / 

i=l '' 




After the proper change of basis, we obtain: 


( 2 ) 

( 3 ) 

( 4 ) 


P{p) 


Nn 


lim A dfj, log Z{ij, — ie) 


( 5 ) 


with: 


z(m) 



d(()i exp (PH (0)) 


( 6 ) 
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^(0) = —fi <p^<p + (f)^J(f) 


This defines a statistics for the fields 0: 


1 r ^ 

= V / n • exp(i?/(0)) 

i=i 


We now exclude the kth field from the statistics and define: 


1 f ^ 

{•)\k = ^ Yld(j)i • exp (mxki^^xk)) 
\k J -/i 


where Z\k, (p\k ci-nd 'H\k are defined like before except the k-th component is removed. 
With this definition, we have: 


{4>l) = ^ / d^kdhk 4>l ( 6{hk - 


The statistics of the fields being gaussian, (6ihk — Jki4>i) ) is gaussian and we 


can write: 


{4>k) = Ykj ^ ^ 


where Zk is the proper normalization constant and: 

<^1 = -2i {hl)\t = -2i ^ ^ (12) 

It can be shown that the second moment involving two different fields is negligible in the 
the large N limit in the expression (one needs to reproduce the same calculation as before 
but excluding fields k and 1 and computing 0^0;). 


We have then: 


4 = 


We can here apply the central limit theorem and obtain: 

N 

al = 

i^k 

At the same time, we can easily derive (0^) as a function of fi and ak- 


k2\ _ d 


{4>k) 


(15) 
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And performing the gaussian integral, we obtain: 


Zh = 


TT 




(16) 


This gives: 


{<A) = 


2/i + all2 


(17) 


By plugging the previously computed we obtain: 


{4>l) = 


2)1-ia^ 


( 18 ) 


Considering now that (pj is self-averaging, we obtain the self-consistent equation: 


IN 


2)1 - ia^(N - l)W))n-i 


( 19 ) 


This relationship ensures that in the limit N —)■ -|-cx): 


VN{{(j)^))N{2fi/^/N - ia^y/NiP?))) + i = 0 


( 20 ) 


((P(/^))) = -limAf((02))iv 

TT 


( 21 ) 


((p(/i))) = 0 if/r>yiVcT 

ip{fj))) = y]|y2\/lVcT2 -/i2 if ^ 


This is the Wigner semi-circle law [1]. 


( 22 ) 


SPECTRUM DENSITY AROUND THE CRITICAL POINT 

In this section, we perturbatively expand the equations (20) of the main text around the 
critical point. To lighten the notation, the averaging brackets shall be implicit. The system 
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of equations writes: 


ki^,n \^? + + mi)^) 

72^Y^ p{k,l,m,n)Jl^ {m-ip}+ n-ii)l) 

ki^,n + + 

,2^sp (l{k,l,m,n)Jl^{k'il)‘f + l'il)l) 
ki^,n \A^ + JLiki^f + + n^l) 

72^Y^ (l{k,l,m,n)J^^{m'ipf + nAl) 

ki^An \z\^ + + nAlAk^f + lipl) 


where V’e = An+b/A f = A'li+i! Jinh- critical value of z, when the spectrum 


density becomes zero, the helds variances tends continuously toward zero. 

order, we obtain; 

At the zero-th 

A‘f= p{.k,l,m,n)J^^^{k^‘l + lijl)/\z\‘^ 

(24a) 



Ae= Y (li^A,m,n)J^^{kij^ + lijl)/\z\‘^ 

(24b) 

k,l,m,n 


Ai= Y Pi^A,m,n)Jl^{mAi +nAl)/\z\‘^ 

(24c) 

k^l,m,n 


It gives the critical value for \z\^ with the held variance ratios: 


Zc ~ Ji h'inh(l f) ~k PexJef 

(25a) 

_ PexJe 

(25b) 

R = ^ zz 

Jf 

(25c) 


(25d) 


(23a) 

(23b) 

(23c) 

(23d) 
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At the first order, one can obtain the spectrum density close to the critical value of z. We 
write the variables at the hrst order: 


kP = + e 
= P 

= -Ro + Cr = — -^ + Cr 

~ ~ ~ 

— Po P 


With these notations, the self-consistent equations (23) write: 


= eRinh^fnh + (^inh^ + (1 “ + ff^L^^P (26a) 

fl^exJinhf^R = + (^c + (1 ~ + /hex</ex)-P (26b) 

~h'inh'^inh(l ~ f)^R — CRex'^exh'inh"^inh + h'ex<^ex(h'ex^c A (1 — /)hinh'^inh + /h'ex'^ex)-^ (26c) 

~h'inh'^inh(^ ~ f)^R ~ ^Rinh'^inh"^ex A /iex'^ex(^c A (1 — /)h'inh'^inh A fIJ,exJex)P (26d) 


By solving these equations, we obtain: 


P = - 
cr = R 
cr = R 


^cRinh'^inh 


((1 - + ^eJJL){4 + (1 - f)f^inhJth + Z/^ex^x) 

- Rinh^gih) + " f)JL + RL/4x) 

° (R?nh(l - f)Jth + RL/4x)(^c + (1 - f)RmhJth + /Rex^Zx) 

h'ex'^ex h'inh'Znh 


(1 - /)RL4h 


+ /rL-^Zx 


(27a) 

(27b) 

(27c) 


This allows us to hnd the spectrum density at the critical point: 




k^l,m,n 


+q{k,l,m,n)f 1 


e -|- J^^{k -\- lRQ)(^Tn -\- uRq) 


(28) 


P(^) 


-P 


^P'inh'Zji 


inh 


>r((l - f)l^inhJ‘nh + + (1 - /)ftnhJ*h + 


(29) 


( 30 ) 
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The ratio between the dense and the sparse expression writes: 

(1 /)h'inh'7ijih /hex'Tf 


Pden 


P{Z) 


= 1 + 


4 

ex 


This is a correction of order 1/p in the sparsity. 


( 31 ) 


LOCALIZATION FOR DENSE MATRICES 

We plot here the same graph as the one presented in the main article bnt for dense 
matrices: 


FIG. 1. 


IPR 
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PERCOLATION THRESHOLD FOR THE TWO-POPULATIONS NETWORK 

We follow here the formalism adopted in [2]. Gn and Gei denote the generating fnnctions 
of the nnmber of edges going ont of respectively the inhibitory and excitatory nodes. For 
the inhibitory generating fnnction: 

+ 00 

Gio{x) = '^qikX^ (32) 

k=0 

where qik is the probability distribntion of the nnmber of ontgoing synapses of an inhibitory 
neuron. In the main text, we use a Poisson distribution: 


=Pci.h(*:) = lithe ‘“•"/k'. 


(33) 







We now consider the cluster sizes which are the number of nodes (neurons) connected to¬ 
gether by inhibitory or excitatory connections. More specihcally, we choose a connection 
at random in the network and consider its output neuron. We then consider the output 
neurons of this neuron, then in turn their output neurons and so on. The cluster size of the 
initial connection is dehned as the number of neuron reached this way. We denote by the 
probability of a cluster to be of size s and consider the associated generating function: 

+00 

Hi{x) = '^qisx" (34) 

s=0 

If we ignore the possibility of a giant cluster containing many recurrent connections, we can 
assume that the clusters are tree like, in the limit of large N. It is then possible to obtain 
a ”Dyson-like” self-consistent relationship between Hi, Gio and Geo'- 

Hi{x) = x{fqeo + (1 - f)qio)+x{fqei + (1 - f)qii)Hi{x) + (35) 

"V* "V* 

no outgoing edge one outgoing edge 

x{fqe2 + {l-f)qi2)Hi{xf^ - (36) 

"-V-' 

two outgoing edges 

= x{fGeo{Hi{x) + (1 - f)G,o{Hi{x))) (37) 


If we now consider the average number of neurons contained in a cluster, we obtain: 

1 

1-/G;o(1)-(1-/)G'o(1) 

The so-called percolation transition, appears when (s) diverges and a giant cluster appears. 
This transition occurs when (/GgQ(l) -1- (1 — /)G'o(l)) = 1, or when the average number 
of outgoing connections from a neuron, irrespective of its inhibitory of excitatory nature is 
one. 
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